Metabarcoding Singapore
OTUs at 98%
Metabarcoding Singapore
OTUs at 98%
- 1 Aim
- 2 Analysis
- 2.1 Assignment of 98% OTUs based on PR2 database
- 2.2 Phyloseq analysis
- 2.2.1 Create phyloseq files for euk after filtering the data
- 2.2.2 Create phyloseq files for proks
- 2.2.3 Create a list for the auto and hetero phyloseq files
- 2.2.4 Bar plot of divisions per station for auto and hetero
- 2.2.5 Bar plot of class per station for auto and hetero
- 2.2.6 Compare by Straight, Site, Moonsoon
- 2.2.7 Main genera for different division of EUkaryote Autotrophs
- 2.2.8 Main species of Mamiellophyceae and Diatoms (Eukaryotes)
- 2.2.9 Heatmaps
- 2.2.10 NMDS
- 2.2.11 Network analysis
Initialize. This file defines all the necessary libraries and variables
source('Metabarcoding Singapore_init.R', echo=FALSE)1 Aim
- Assign and analyze eukaryotes for Singapore metabarcoding data (98%). Do some analyzes with the prokaryotes too…
2 Analysis
2.1 Assignment of 98% OTUs based on PR2 database
2.1.1 Read the data
# Read the sample table
sample_table <- read_excel("../qiime/Singapore otu_table_w_tax_0.98.xlsx", sheet = "samples")
# Read the otu table
otu_table <- read_excel("../qiime/Singapore otu_table_w_tax_0.98.xlsx", sheet = "otu_table")
# Clean up the taxonomy
otu_table <- otu_table %>% mutate(taxo_clean = str_replace_all(taxonomy, "D_[0-9]+__",
"")) %>% separate(col = taxo_clean, into = str_c("taxo", c(1:7)), sep = "; ")
# Read the sequences
otu_sequences <- readAAStringSet("../qiime/otu_rep_98.fna")
otu_sequences.df <- data.frame(names = names(otu_sequences), sequence = as.character(otu_sequences))
# Split the otu sequences name into the otu id given by qiime and the name
# of the representative sequence Remove the primers
fwd_length = 20
rev_length = 15
otu_sequences.df <- otu_sequences.df %>% separate(col = names, into = c("otu_id_qiime",
"otu_rep_seq"), sep = " ") %>% mutate(sequence = str_sub(sequence, start = fwd_length +
1, end = -rev_length - 1))
otu_table <- left_join(otu_table, otu_sequences.df)
# Write a fasta file for blast with all taxonomy roups
otu_sequences <- otu_table %>% transmute(sequence = sequence, seq_name = otu_id)
fasta_write(otu_sequences, file_name = "../qiime/otu_rep_98_all.fasta", compress = FALSE,
taxo_include = FALSE)[1] TRUE
2.1.2 Only keep the eukaryotes in the OTU file
otu_table_euk <- otu_table %>% filter(str_detect(taxonomy, "Eukaryota|Unassigned"))
# Write the fasta file file
otu_sequences_euk <- otu_table_euk %>% transmute(sequence = sequence, seq_name = otu_id)
fasta_write(otu_sequences_euk, file_name = "../qiime/otu_rep_98_euk.fasta",
compress = FALSE, taxo_include = FALSE)[1] TRUE
2.1.3 Use dada2 to reassign to PR2
dada2_assign(seq_file_name = "../qiime/otu_rep_98_euk.fasta")2.1.4 Read the PR2 assignement and merge with initial otu table
otu_euk_pr2 <- read_tsv("../qiime/otu_rep_98_euk.dada2.taxo")
otu_euk_pr2_boot <- read_tsv("../qiime/otu_rep_98_euk.dada2.boot") %>% rename_all(funs(str_c(.,
"_boot"))) %>% rename(seq_name = seq_name_boot)
otu_euk_pr2 <- left_join(otu_euk_pr2, otu_euk_pr2_boot) %>% rename(otu_id = seq_name)
otu_table_final <- left_join(otu_table, otu_euk_pr2) %>% select(otu_id, otu_id_qiime,
taxo1:taxo7, kingdom:species_boot, nseq, STJ11.S23:PR06.S10, otu_rep_seq,
sequence)
write_tsv(otu_table_final, "../qiime/otu_table_98_final.tsv")2.2 Phyloseq analysis
2.2.1 Create phyloseq files for euk after filtering the data
Filter the euk data to remove the low bootstraps values (threshold : bootstrap > 90% at the supergroup level) and create a phyloseq file
Note the bootstrap threshold had to be higher for 98% compared to 97% (90% vs 65%)
otu_table_euk_final <- otu_table_final %>% filter(supergroup_boot > 90)
str_c("Number of euk otus : ", nrow(otu_table_euk_final))[1] "Number of euk otus : 587"
otu_mat <- otu_table_euk_final %>% select(otu = otu_id, STJ11.S23:PR06.S10)
tax_mat <- otu_table_euk_final %>% select(otu = otu_id, kingdom:species)
samples_df <- sample_table %>% rename(sample = sample_id)
row.names(otu_mat) <- otu_mat$otu
otu_mat <- otu_mat %>% select(-otu)
row.names(tax_mat) <- tax_mat$otu
tax_mat <- tax_mat %>% select(-otu)
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select(-sample)
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
singa_euk <- phyloseq(OTU, TAX, samples)
singa_eukphyloseq-class experiment-level object
otu_table() OTU Table: [ 587 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 587 taxa by 8 taxonomic ranks ]
Break up into photosynthetic and non-photosynthetic (Metazoa are removed)
singa_photo <- subset_taxa(singa_euk, (division %in% c("Chlorophyta", "Dinoflagellata",
"Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta", "Cercozoa")) &
!(class %in% c("Syndiniales", "Sarcomonadea")))
singa_hetero <- subset_taxa(singa_euk, !(division %in% c("Chlorophyta", "Dinoflagellata",
"Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta", "Cercozoa", "Metazoa")) |
(class %in% c("Syndiniales", "Sarcomonadea")))
sprintf("Euk Photo - all")[1] "Euk Photo - all"
singa_photophyloseq-class experiment-level object
otu_table() OTU Table: [ 178 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 178 taxa by 8 taxonomic ranks ]
sprintf("Euk Hetero - all")[1] "Euk Hetero - all"
singa_heterophyloseq-class experiment-level object
otu_table() OTU Table: [ 260 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 260 taxa by 8 taxonomic ranks ]
Normalize number of reads in each sample using median sequencing depth.
total_photo = median(sample_sums(singa_photo))
sprintf("The number of reads used for normalization of autotrophs is %.0f",
total_photo)[1] "The number of reads used for normalization of autotrophs is 383"
standf = function(x, t = total_photo) (t * (x/sum(x)))
singa_photo = transform_sample_counts(singa_photo, standf)
total_hetero = median(sample_sums(singa_hetero))
sprintf("The number of reads used for normalization of heterotrophs is %.0f",
total_hetero)[1] "The number of reads used for normalization of heterotrophs is 100"
standf = function(x, t = total_hetero) (t * (x/sum(x)))
singa_hetero = transform_sample_counts(singa_hetero, standf)Remove taxa which are in low abundace (< 0.10 in any given sample) and normalize again…
contrib_min <- 0.1
singa_photo_abundant <- filter_taxa(singa_photo, function(x) sum(x > total_photo *
contrib_min) > 0, TRUE)
singa_hetero_abundant <- filter_taxa(singa_hetero, function(x) sum(x > total_hetero *
contrib_min) > 0, TRUE)
total = median(sample_sums(singa_photo_abundant))
sprintf("The number of reads used for normalization of abundant autotrophs is %.0f",
total)[1] "The number of reads used for normalization of abundant autotrophs is 290"
standf = function(x, t = total) (t * (x/sum(x)))
singa_photo_abundant = transform_sample_counts(singa_photo_abundant, standf)
total = median(sample_sums(singa_hetero_abundant))
sprintf("The number of reads used for normalization of abundant heterotrophs is %.0f",
total)[1] "The number of reads used for normalization of abundant heterotrophs is 52"
standf = function(x, t = total) (t * (x/sum(x)))
singa_hetero_abundant = transform_sample_counts(singa_hetero_abundant, standf)
sprintf("Euk Photo - abundant")[1] "Euk Photo - abundant"
singa_photo_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 26 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 26 taxa by 8 taxonomic ranks ]
sprintf("Euk Hetero - abundant")[1] "Euk Hetero - abundant"
singa_hetero_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 46 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 46 taxa by 8 taxonomic ranks ]
2.2.2 Create phyloseq files for proks
otu_table_prok <- otu_table %>% filter(taxo1 %in% c("Bacteria", "Archaea"))
otu_mat <- otu_table_prok %>% select(otu = otu_id, STJ11.S23:PR06.S10)
tax_mat <- otu_table_prok %>% select(otu = otu_id, taxo1:taxo7) %>% rename(kingdom = taxo1,
supergroup = taxo2, division = taxo3, class = taxo4, order = taxo5, family = taxo6,
genus = taxo7) %>% mutate(species = NA)
samples_df <- sample_table %>% rename(sample = sample_id)
row.names(otu_mat) <- otu_mat$otu
otu_mat <- otu_mat %>% select(-otu)
row.names(tax_mat) <- tax_mat$otu
tax_mat <- tax_mat %>% select(-otu)
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select(-sample)
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
singa_prok <- phyloseq(OTU, TAX, samples)
sprintf("Prok - all")[1] "Prok - all"
singa_prokphyloseq-class experiment-level object
otu_table() OTU Table: [ 1476 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 1476 taxa by 8 taxonomic ranks ]
Normalize number of reads in each sample using median sequencing depth.
total_prok = median(sample_sums(singa_prok))
sprintf("The number of reads used for normalization of prokaryotes is %.0f",
total_prok)[1] "The number of reads used for normalization of prokaryotes is 7334"
standf = function(x, t = total_prok) (t * (x/sum(x)))
singa_prok = transform_sample_counts(singa_prok, standf)Remove otus which are in low abundace (< 0.05 in any given sample) and normalize again
contrib_min <- 0.05
singa_prok_abundant <- filter_taxa(singa_prok, function(x) sum(x > total_prok *
contrib_min) > 0, TRUE)
total = median(sample_sums(singa_prok_abundant))
sprintf("The number of reads used for normalization of abundant prokaryotes is %.0f",
total)[1] "The number of reads used for normalization of abundant prokaryotes is 3187"
standf = function(x, t = total) (t * (x/sum(x)))
singa_prok_abundant = transform_sample_counts(singa_prok_abundant, standf)
sprintf("Prok - abundant")[1] "Prok - abundant"
singa_prok_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 42 taxa and 64 samples ]
sample_data() Sample Data: [ 64 samples by 5 sample variables ]
tax_table() Taxonomy Table: [ 42 taxa by 8 taxonomic ranks ]
2.2.3 Create a list for the auto and hetero phyloseq files
ps_list <- list(ps = c(singa_photo, singa_hetero, singa_prok), title = c("Eukaryotes - Autotrophs - all OTUs",
"Eukaryotes - Heterotrophs - all OTUs", "Prokaryotes - all OTUs"))
ps_list_abundant <- list(ps = c(singa_photo_abundant, singa_hetero_abundant,
singa_prok_abundant), title = c("Eukaryotes - Autotrophs - abundant OTUs (> 10%)",
"Eukaryotes - Heterotrophs - abundant OTUs (> 10%)", "Prokaryotes - abundant OTUs (> 5%)"))2.2.4 Bar plot of divisions per station for auto and hetero
for (i in 1:3) {
p <- plot_bar(ps_list$ps[[i]], fill = "division") + geom_bar(aes(color = division,
fill = division), stat = "identity", position = "stack") + ggtitle(str_c("Division level - ",
ps_list$title[[i]]))
print(p)
}2.2.5 Bar plot of class per station for auto and hetero
Only consider the abundant taxa
for (i in 1:3) {
p <- plot_bar(ps_list_abundant$ps[[i]], fill = "class") + geom_bar(aes(color = class,
fill = class), stat = "identity", position = "stack") + ggtitle(str_c("Class level - ",
ps_list_abundant$title[[i]]))
print(p)
}2.2.6 Compare by Straight, Site, Moonsoon
for (factor in c("strait", "site", "moonsoon")) {
for (i in 1:3) {
ps_aggregate <- merge_samples(ps_list_abundant$ps[[i]], factor)
p <- plot_bar(ps_aggregate, fill = "division") + geom_bar(aes(color = division,
fill = division), stat = "identity", position = "stack") + ggtitle(str_c(ps_list_abundant$title[[i]],
" - ", factor))
print(p)
}
}2.2.7 Main genera for different division of EUkaryote Autotrophs
for (one_division in c("Chlorophyta", "Dinoflagellata", "Ochrophyta")) {
ps_subset <- subset_taxa(singa_photo_abundant, division %in% one_division)
p <- plot_bar(ps_subset, x = "genus", fill = "genus", facet_grid = strait ~
moonsoon) + geom_bar(aes(color = genus, fill = genus), stat = "identity",
position = "stack") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5,
hjust = 1)) + ggtitle(str_c(one_division, " - Abundant OTUs"))
print(p)
}2.2.8 Main species of Mamiellophyceae and Diatoms (Eukaryotes)
for (one_class in c("Mamiellophyceae", "Bacillariophyta")) {
ps_subset <- subset_taxa(singa_photo_abundant, class %in% one_class)
p <- plot_bar(ps_subset, x = "species", fill = "species", facet_grid = strait ~
moonsoon) + geom_bar(aes(color = species, fill = species), stat = "identity",
position = "stack") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5,
hjust = 1)) + ggtitle(str_c(one_class, " - Abundant OTUs"))
print(p)
}2.2.9 Heatmaps
for (i in 1:2) {
p <- plot_heatmap(ps_list_abundant$ps[[i]], method = "NMDS", distance = "bray",
taxa.label = "species", taxa.order = "division", low = "beige", high = "red",
na.value = "beige", title = ps_list_abundant$title[[i]])
print(p)
}for (i in c(3)) {
p <- plot_heatmap(ps_list_abundant$ps[[i]], method = "NMDS", distance = "bray",
taxa.label = "family", taxa.order = "division", low = "beige", high = "red",
na.value = "beige", title = ps_list_abundant$title[[i]])
print(p)
}2.2.10 NMDS
All OTUs
for (i in 1:3) {
singa.ord <- ordinate(ps_list$ps[[i]], "NMDS", "bray")
print(singa.ord)
p <- plot_ordination(ps_list$ps[[i]], singa.ord, type = "samples", color = "moonsoon",
shape = "strait", title = ps_list$title[[i]]) + geom_point(size = 3)
print(p)
}Square root transformation
Wisconsin double standardization
Run 0 stress 0.14
Run 1 stress 0.15
Run 2 stress 0.15
Run 3 stress 0.17
Run 4 stress 0.16
Run 5 stress 0.15
Run 6 stress 0.16
Run 7 stress 0.15
Run 8 stress 0.16
Run 9 stress 0.15
Run 10 stress 0.16
Run 11 stress 0.15
Run 12 stress 0.16
Run 13 stress 0.16
Run 14 stress 0.14
Run 15 stress 0.16
Run 16 stress 0.16
Run 17 stress 0.16
Run 18 stress 0.16
Run 19 stress 0.15
Run 20 stress 0.15
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.14
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.18
Run 1 stress 0.19
Run 2 stress 0.18
Run 3 stress 0.18
Run 4 stress 0.17
... New best solution
... Procrustes: rmse 0.029 max resid 0.16
Run 5 stress 0.18
Run 6 stress 0.19
Run 7 stress 0.18
Run 8 stress 0.18
Run 9 stress 0.18
Run 10 stress 0.19
Run 11 stress 0.2
Run 12 stress 0.18
Run 13 stress 0.19
Run 14 stress 0.18
Run 15 stress 0.18
Run 16 stress 0.18
Run 17 stress 0.18
Run 18 stress 0.18
Run 19 stress 0.18
Run 20 stress 0.18
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.17
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.089
Run 1 stress 0.11
Run 2 stress 0.1
Run 3 stress 0.095
Run 4 stress 0.11
Run 5 stress 0.1
Run 6 stress 0.09
Run 7 stress 0.095
Run 8 stress 0.094
Run 9 stress 0.11
Run 10 stress 0.1
Run 11 stress 0.098
Run 12 stress 0.11
Run 13 stress 0.1
Run 14 stress 0.11
Run 15 stress 0.11
Run 16 stress 0.1
Run 17 stress 0.11
Run 18 stress 0.11
Run 19 stress 0.1
Run 20 stress 0.09
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.089
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Abundant OTUs
for (i in 1:3) {
singa.ord <- ordinate(ps_list_abundant$ps[[i]], "NMDS", "bray")
print(singa.ord)
p <- plot_ordination(ps_list$ps[[i]], singa.ord, type = "samples", color = "moonsoon",
shape = "strait", title = ps_list_abundant$title[[i]]) + geom_point(size = 3)
print(p)
p <- plot_ordination(ps_list$ps[[i]], singa.ord, type = "taxa", color = "division",
shape = "division", title = ps_list_abundant$title[[i]]) + geom_point(size = 3)
print(p)
}Square root transformation
Wisconsin double standardization
Run 0 stress 0.14
Run 1 stress 0.18
Run 2 stress 0.16
Run 3 stress 0.18
Run 4 stress 0.17
Run 5 stress 0.18
Run 6 stress 0.17
Run 7 stress 0.18
Run 8 stress 0.18
Run 9 stress 0.18
Run 10 stress 0.17
Run 11 stress 0.16
Run 12 stress 0.18
Run 13 stress 0.17
Run 14 stress 0.18
Run 15 stress 0.17
Run 16 stress 0.16
Run 17 stress 0.15
Run 18 stress 0.17
Run 19 stress 0.17
Run 20 stress 0.16
*** No convergence -- monoMDS stopping criteria:
18: stress ratio > sratmax
2: scale factor of the gradient < sfgrmin
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.14
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Wisconsin double standardization
Run 0 stress 0.24
Run 1 stress 0.25
Run 2 stress 0.25
Run 3 stress 0.25
Run 4 stress 0.26
Run 5 stress 0.25
Run 6 stress 0.26
Run 7 stress 0.26
Run 8 stress 0.26
Run 9 stress 0.25
Run 10 stress 0.26
Run 11 stress 0.25
Run 12 stress 0.25
Run 13 stress 0.26
Run 14 stress 0.26
Run 15 stress 0.27
Run 16 stress 0.24
... Procrustes: rmse 0.012 max resid 0.06
Run 17 stress 0.25
Run 18 stress 0.24
Run 19 stress 0.25
Run 20 stress 0.27
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(veganifyOTU(physeq))
Distance: bray
Dimensions: 2
Stress: 0.24
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(veganifyOTU(physeq))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.12
Run 1 stress 0.12
Run 2 stress 0.13
Run 3 stress 0.14
Run 4 stress 0.13
Run 5 stress 0.12
Run 6 stress 0.12
Run 7 stress 0.13
Run 8 stress 0.13
Run 9 stress 0.14
Run 10 stress 0.13
Run 11 stress 0.12
Run 12 stress 0.13
Run 13 stress 0.13
Run 14 stress 0.12
... New best solution
... Procrustes: rmse 0.0024 max resid 0.018
Run 15 stress 0.14
Run 16 stress 0.13
Run 17 stress 0.12
Run 18 stress 0.12
Run 19 stress 0.13
Run 20 stress 0.14
*** No convergence -- monoMDS stopping criteria:
18: stress ratio > sratmax
2: scale factor of the gradient < sfgrmin
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.12
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
2.2.11 Network analysis
for (i in 1:2) {
p <- plot_net(ps_list_abundant$ps[[i]], distance = "(A+B-2*J)/(A+B)", type = "taxa",
maxdist = 0.4, color = "class", point_label = "genus") + ggtitle(ps_list_abundant$title[[i]])
print(p)
}for (i in c(3)) {
p <- plot_net(ps_list_abundant$ps[[i]], distance = "(A+B-2*J)/(A+B)", type = "taxa",
maxdist = 0.4, color = "class", point_label = "family") + ggtitle(ps_list_abundant$title[[i]])
print(p)
}